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Abstract 

(N 

^ I The onset of synchronization in networks of networks is investigated. SpecificaUy, we consider 

'NT 

, frequencies are drawn from a given distribution, and each population has its own such distribution. 

o 



networks of interacting phase oscihators in which the set of osciUators is composed of several dis- 
tinct populations. The oscillators in a given population are heterogeneous in that their natural 



The coupling among the oscillators is global, however, we permit the coupling strengths between 
ij - the members of different populations to be separately specified. We determine the critical condi- 

. tion for the onset of coherent collective behavior, and develop the illustrative case in which the 

oscillator frequencies are drawn from a set of (possibly different) Cauchy-Lorentz distributions. 
One motivation is drawn from neurobiology, in which the collective dynamics of several interacting 
populations of oscillators (such as excitatory and inhibitory neurons and glia) are of interest. 
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In recent years, there has been considerable interest in networks of interacting systems. 
Researchers have found that an appropriate description of such systems involves an under- 
standing of both the dynamics of the individual oscillators and the connection topology of 
the network. Investigators studying the latter have found that many complex networks have 
a modular structure involving motifs communities , layers [^], or clusters [5|. For 
example, recent work has shown that as many kinds of networks (including isotropic homo- 
geneous networks and a class of scale-free networks) transition to full synchronization, they 
pass through epochs in which well-defined synchronized communities appear and interact 
with one another Q| . Knowledge of this structure, ^nd the dynamical behavior it supports, 
informs our understanding of biological |6i], social [7], and technological networks 8|]. 

Here we consider the onset of coherent collective behavior in similarly structured systems 
for which the dynamics of the individual oscillators is simple. In seminal work, Kuramoto 
analyzed a mathematical model that illuminates the mechanisms by which synchronization 
arises in a large set of globally-coupled phase oscillators [9]. An important feature of Ku- 
ramoto's model is that the oscillators are heterogeneous in their frequencies. And, although 
these mathematical results assume global coupling, they have been fruitfully applied to fur- 
ther our understanding of systems of fireflies, arrays of Josephson junctions, electrochemical 
oscillators, and many other cases 10 1. 

In this work, we study systems of several interacting Kuramoto systems, i.e., networks 
of interacting populations of phase oscillators. Our motivation is drawn not only from the 
applications listed above (e.g., an amusing application might be interacting populations of 
fireflies, where each population inhabits a separate tree), but also from other biological sys- 
tems. Rhythms arising from coupled cell populations are seen in many of the body's organs 
(including the heart, the pancreas, and the kidneys, to name but a few), all of which interact. 
For example, the circadian rhythm influences many of these systems. We draw additional 
motivation from consideration of how neuronal tissue is organized. Although we do not 
consider neuronal systems specifically in this paper, we note that heterogeneous ensembles 
of neurons often exhibit a "network-of-networks" topology. At the cellular level, populations 
of particular kinds of neurons (e.g., excitatory neurons) interact not only among themselves, 
but also with populations of other distinct neuronal types (e.g., inhibitory neurons). At 



a higher level o 
another as well 



organization, various anatomical regions of the brain interact with one 
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Although our network is simple, it is novel in that we include heterogeneity at several 
levels. Each population consists of a collection of phase oscillators whose natural frequencies 
are drawn at random from a given distribution. To allow for heterogenity at the population 
level, we let each population have its own such frequency distribution. In addition, our 
system is heterogeneous at the coupling level as well: we consider global coupling such 
that the coupling strengths between the members of different populations can be separately 
specified. 

The assumpion of global (but population-weighted) coupling permits an analytic deter- 
mination of the critical condition for the onset of coherent collective behavior, as we will 
show. While this assumption may not strictly apply in some applications, our results provide 
insight into the behavior of networks of similar topology even if the connectivity is less than 
global. 

We begin by specifying our model and deriving our main results. We then discuss several 
illustrative examples. 

Consider first a two-species Kuramoto model. We label the separate species 9 and (p and 
assume that there are No and A^^ such oscillators in each population, respectively. Thus, 
the system equations are 

dO- k k ^'^ 

—^ = ujei + sin{9j - 9i - aee) + ^2 ■^^^(^i - - "e^), 

-—^ = uj^i + 7]^ V siniOj a^g) + r/-^ V sin{(j)j -(pi- a^^). 

Here, rj is an overall coupling strength, the a's provide additional heterogeneity in the 
coupling functions, and the matrix 

K = I V I 

kd^ft kA, 



defines the connectivity among the populations 

For arbitrarily many different species, let a range over the various population symbols 
6,(j), . . . with the understanding that depending on the context, a may represent either a 
label (when subscripted) or a variable. Thus, we have 



The ujai are the constant natural frequencies of the oscillators when uncoupled, and are 



distributed according to a set of time- independent distribution functions Ga{oJa) [12|. 
We define the usual Kuramoto order parameter for each population, i.e., 



Here, describes the degree of synchronization in each population, and ranges from to 1. 
Using this, the above equations can be expressed 



dt 



uJai + ^r]kaa'ra'Sin{tp„> - <Ji - a„^>). 



Assuming that the are very large, we solve for the onset of coherent collective be- 
havior by using a distribution function approach. Thus, instead of discrete indicies, we 
imagine continua of oscillators described by distribution functions Fo-(cr, cUo-, t) such that 
Ffj{(T,uJr7,t)daduJcr is the fraction of a-oscillators whose phases and natural frequencies lie in 
the infinitesimal volume element daduj^ at time t. Note that in the N^j oo limit, 

Ga{i^a) = / F^{a,uj„,t)da, 
Jo 

and the order parameters are given by 

r^e'^- = [ [ F^e'^dadu^. (2) 



In this context, the distribution functions satisfy the equations of continuity, i.e. 

dFfj -t , ^ da ^. 
_ + V . (F,-a) = 0, 

and if we write F(j{cr,u^,t) = f(j{a,uJa,t)Ga{uj^), we have 



dt da 



/. = 0- (3) 



The incoherent state has a distributed uniformly over [0, 27r], so that fa = 1/2tt and = 
0. These satisfy Eq. trivially. We test the stability of this solution by perturbing it. Note 
that a perturbation to leads to a perturbation of r^, and we expect these perturbations 
to either grow or shrink exponentially in time, depending on the overall coupling strength r]. 
The marginally stable case occurs at a critical value rj^, at which coherent collective behavior 
emerges. 



Thus, we write fa = l/27r + (5/o-)e** and = (5ro-)e^*, where (Sfa) and (Sr^) are smalL 
Inserting the first of these into the continuity equation (Eq. ([3])), and keeping only first-order 



terms, we obtain 



da 

The solution to this equation is 

1 



d 1 



a' 



+ 



S — lUJn 



s + iu„ 



(4) 



Consistency demands that the perturbations {5 fa) and {5ra) be related to each other 
via the order parameter equation, Eq. ([2]). This yields our main result, as follows. Eq. ([2]) 
becomes 



st ilpa 



Ga(i^a) 



2tt 



27V 



{6fa)e'' e'^daduja. 



The integral involving l/27r is zero. Inserting the solution for {6 fa) from Eq. (jl]), one obtains 



{5ra)e 



1 r 



2 ./_oo S - tUJc 



duJa 



Define the bracketed expression to be ga{s), and define Za = {6ra)e'^'^'^ . Then, we have 

Za = 9a{s)^rikaa'e~"^-''' Za' . 
a' 

Now define the complex quantity kaa' = vkaa'^'^"''"'' ■ Using the Kronecker delta Saa'i the 
above equation can be written 



9 As] 



Za' = 0. 



In matrix notation for the case of two populations labeled 6 and 0, this is 



geis) ■"'"P 



0. 



(5) 



9<#.(s) 



This equation has a non-trivial solution if the determinant of the matrix is zero. This 
condition determines the growth rate s in terms of ?7, kaa') and the parameters of the 
distributions Gai^Ja) (via ga{s)). 
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To illustrate the resulting behavior, we note that gcr{s) can be evaluated in closed form 
for a Cauchy-Lorentz distribution 

where Q^, the mean frequency of population a, and the half-width-at-half-maximum A^- are 
both real, and Ao- is positive. One obtains 

1 



2{s + A^ + in^). 

Using this, the determinant condition for the two-population case reduces to 

[kee - 2(s + Ag + iQg)] [k^^ - 2(s + + iQ^)] - kg^k^g = 0. 

For simplicity, we set the phase angles ao-o-' to zero for the remainder of this paper [l^ . 
In this case, the matrix elements fco-o-' are purely real, so that k„„/ = rjk^^r . The determinant 
condition then becomes 

[rikgg - 2(s + Ag + iQg)] [r]k^^ - 2{s + A,^ + iQ^)] - rfkg^k^g = 0. (7) 

Notice that if r/ = 0, then s = —A^ — iQ^, indicating that the incoherent state is stable 
for zero coupling (since — A^- is negative). We imagine increasing (or decreasing) rj until s 
crosses the imaginary axis at a critical value r^*. At this point, the incoherent state loses 
stability and coherent collective behavior emerges in the ensemble. Thus, the critical value rj^, 
can be determined from the determinant condition by setting s = iv, where v is real (so that 
the perturbations are marginally stable), and equating the real and imaginary parts of the 
left side of Eq. ([7]) to zero. This results in two equations which can be solved simultaneously 
for the two (real) unknowns f] and v. 

For our first example, we choose two identical populations; we set Ag = A^ = A and 
Qg = Qfl^ = Q. (A more generic example will follow.) Denoting D = det(K) and T = tr(K), 
we separate the real and imaginary parts of Eq. ([7]) to obtain 

Dt]^ - 2ATr] + 4A2 - A{v + = 
{v + n){AA-riT) = 0. 

One solution to these equations is 

'T±y/T^-AD 



—Q, ?7* = A 



D 



which is vahd for > 4D, since t] is assumed to be real. Notice that the appropriate 
solution as ^ (using the negative sign for T > and the positive sign for T < 0) is 

-_o 

as can be verified by setting D = in Eq. ([8]) directly. Another solution is 

Finally, setting T = in Eq. ([8]) yields = —Vl and 77* = ^^=p for D < 0, and no solution 
for D >0. These results are summarized in Table [H 

Thus, the critical values rj^, are determined by T, D, and A. To illustrate this result, we 
begin by discussing a particular example. Consider the matrix 



K 



1 -1 
1 



which has trace T = 1 and determinant D = 1. This corresponds to case two in Table [H 
from which we find that r/* = 4A/T. Fig. [T] shows the results of a numerical simulation of 
two populations of 10, 000 oscillators each, using A = 1. The order parameters r^r versus r) 
are shown, and we can see that the oscillator populations are incoherent for rj values below 
the predicted critical value r/* = 4, and that they grow increasingly synchronized as 1] is 
increased beyond 77*. 

Next, we examined eight different connectivity matricies K that were chosen to sample 
the various regions in T — D space. Table shows these matricies and the corresponding 
value(s) of r^*. These predictions were tested by numerically calculating the order param- 



eters r„ for a range of coupling values t] ISj]. Results are shown in Fig. [2J The various 



cases are located by letter in the D — T plane according to the trace and determinant of the 
matricies, and the corresponding inset shows the numerically-calculated order parameters 
plotted versus 77, with the predicted critical coupling indicated by a vertical line at 77 = 77^,. 
For example, the inset to case A (with {T,D) = (1, 1)) corresponds to Fig. [T], which was 
discussed above. It can be seen that for all cases, the onset of synchronization occurs as 
predicted. 

Note that more than one prediction for 77* may be specified by our analysis (see Table 
[T]). The solutions closest to 7/ = are the relevant ones, because we expect the incoherent 
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state to lose stability once the first 77* solution is encountered. There are two possible cases 
depending on the sign of D. First, if the two solutions have the same sign, then there is 
only one critical value rj^, (which may be positive or negative depending on the sign of the 
trace) for the onset of synchronization. This occurs for D > Q and T 7^ 0, as can be seen 
in Fig. [21 (Interestingly, for D > and T = 0, synchronization does not occur for any rj.) 
The other case occurs for D < 0, for which the two rj^ solutions have opposite signs. In this 
case, there are two critical values r/* for the onset of synchronization - one on either side of 
7] = 0. This can also be observed in Fig. [2l 

In the more general case in which the various populations have different natural frequency 
distributions, it is not typically possible to describe the onset of synchronization in terms 
of the determinant and trace of the coupling matrix K = k^^i alone. We now consider 
this situation, but retain the Cauchy-Lorentz form of the natural frequency distributions 
for convenience. We manipulate Eq. ([7j) as follows. Let s = iv (i.e., purely imaginary, to 
consider the marginally stable case) and define a = Ag + i{v + Qg) and b = + i{v + fi<^). 
We obtain 

Dr]^ - 2ri{bkgg + ak^^) + Aab = 0. (9) 

There are two unknowns in this equation: rj and v. Eq. ([H]) is a quadratic equation in rj with 
complex coefficients, and we can easily obtain two complex solutions r/*^^'^^ as functions of v. 
Since the critical values rji^'^'' must be real, we solve for the roots of 1 711(1]^^''^^), and evaluate 
^g(^i^(i,2)^ at these roots. This yields the possible critical values r]i^''^\ Typically, these steps 
must be performed symbolically and/or numerically; we used MATLAB® [16]. As before, 
the values of rjl that are closest to zero (on either side) are the relevant values. 

To illustrate, we choose two populations with Cauchy-Lorentz natural frequency distri- 
butions (Eq. (j6l)) with parameters Ag = 1, Qg = 2, = 0.5, and fl^p = 4. We consider 
the same K matricies as before, i.e., those listed in the second column of Table HTl Case E 
is straightforward; the analysis is illustrated in Fig. ([3]) and the numerical verification is in 
Fig. (jl]). Note that since Eq. (Q has complex coefficients, obtaining rj typically involves 
taking the square root of a complex number. Therefore, one must be mindful of branch cuts 
when obtaining symbolic and/or numerical solutions. This is important in the analysis for 
case A, shown in Figs. ([5]) and (EI). Finally, case H, which exhibits no synchronization for 
identical populations for any value of rj, does show synchronization in the present case. The 
analysis is shown in Figs. ([3IH1)- 
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We close by giving an example with three different populations. We choose the same 
Cauchy-Lorentz distributions as above, and add a third with Ap = 1/3 and Qp = 1. We use 
the following K matrix: 

/-II 1 \ 

1-11 

\ 1 1 -ij 

The procedure for deriving 77^, proceeds as above, except that Eq. ([7]) is replaced by a third- 
degree polynomial in rj. Fig. ([9]) shows the imaginary and real parts of the three r] solutions. 
(Note that the branch cuts are more complicated.) The predicted onset of synchronization 
was verified, as shown in Fig. ffTOl) . 

In conclusion, we have described how to determine the onset of coherent collective be- 
havior in systems of interacting Kuramoto systems, i.e., systems of interacting populations 
of phase oscillators with both node and coupling heterogeneity. EB was supported by NIH 
grant R01-MH79502; EO was supported by ONR (Physics) and by NSF grant PHY045624. 



[1] R. Milo et al., Science 298, 824-827 (2002). 

[2] M.E.J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004). 

[3] A. Arenas, A. Di'az-Guilera, and C.J. Perez- Vicente, Phys. Rev. Lett. 96, 114102 (2006); A. 

Arenas, A. Diaz-Guilera, and C.J. Perez- Vicente, Physica D 224, 27-34 (2006). 
[4] M. Kurant and P. Thiran, Phys. Rev. Lett. 96, 138701 (2006). 
[5] X. Wang, L. Huang, Y-C Lai, and C.H. Lai, Phys. Rev. E 76, 056113 (2007). 
[6] C. Zhou, L. Zemanova, G. Zamora, C.C. Hilgetag, and J. Kurths, Phys. Rev. Lett. 97, 238103 

(2006). 

[7] M. Girvan and M.E.J. Newman, Proc. Natl. Acad. Sci. 99 7821-7826 (2002). 
[8] R. Milo et al.. Science 303, 1538-1542 (2004). 

[9] Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, 
edited by H. Araki, Lecture Notes in Physics, Vol. 39 (Springer, Berlin, 1975); Chemical 
Oscillators, Waves and Turbulence (Springer, Berlin, 1984). 
[10] A. T. Winfree, The Geometry of Biological Time (Springer, New York, 1980); S.H. Strogatz, 
Physica D 143, 1 (2000); J.A. Acebron et al.. Rev. Mod. Phys. 77, 137-185. K. Wiesenfeld 



9 



and J.W. Swift, Phys. Rev. E 51, 1020-1025 (1995). I.Z. Kiss, Y. Zhai, and J.L. Hudson, 

Science 296 1676 (2005). 
[11] A similar system of two asymmetrically interacting populations with a particular form of 

coupling matrix K was considered in E. Montbrio, J. Kurths, and B. Blasius, Phys. Rev. E 

70, 056125 (2004). The formulation in the current paper is more general in two important 

aspects: K is arbitrary, and we allow for any number of interacting populations. 
[12] Our system is similar to that studied in J.G. Restrepo, E. Ott, and J.G. Restrepo, Chaos 16, 

015107 (2005). However, our formulation permits different natural frequency distributions for 

each population. 

[13] The bracketed expression is valid for Re{s) > and may be analytically continued into the 
region where Re{s) < 0. See E. Ott, P. So, E. Barreto, and T. Antonsen, Physica D 173, 
226-258 (2002); E. Ott, Chaos in Dynamical Systems, second edition, Cambridge University 
Press, 2002, p. 240. 

[14] Many interesting states require Oa-a' 7^ 0, such as the chimera state observed in D.M. Abrams 
and S.H. Strogatz, Int. J. Bif. Chaos 16 #1, 21-37 (2006). 

[15] Simulations used fourth-order Runge-Kutta with a timestep of 0.01 seconds, = 10, 000 or 
50, 000, and parameters as noted. The system was initialized in the incoherent state and an 
initial transient was discarded. The order parameters were then averaged over the subsequent 
10 seconds. Because the standard deviations over this interval were small, no error bars were 
plotted. 

[16] MATLAB®, The MathWorks, Inc., Natick, MA (jhttp : //www . mathworks . com/p . 



Figures 



10 




FIG. 1: Numerical calculation of the order parameters (•, A) versus r] for case A in Table Ull The 
vertical line corresponds to the predicted value r/^, = 4. The data point nearest r]^: is at i] = 4.15. 
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FIG. 2: Numerical simulations using the matricies listed in Table (jll]) for identical populations. 
The letters indicate the placement of each case in the T — D plane, and the corresponding insets 
show numerical calculations of the order parameters (•, A) versus r] for that case (in all cases, 
r/ = is in the center of the horizontal axis). Vertical lines in the insets indicate the predicted 
value(s) r/* listed in Table (jTl|) for the onset of coherent collective behavior. In all cases, we used 
A = 1. Note that for D > 0, there is one value of r]^^ whose sign corresponds to the sign of the 
trace T. If Z) > and T = 0, then synchronization is not possible for any coupling strength. For 
Z) < 0, there are two values of r/*: one positive, and one negative. 
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-4 -3 -2 -1 -4 -3 -2 -1 

V V 

FIG. 3: Case E, different populations. Tlie left panel shows Im{ri^^''^^), with the vertical lines 
identifying roots at vi = —4.024 and V2 = —1.722. The right panel shows Re{r]^^''^^); values at the 
roots found above are indicated by horizontal lines, yielding rji^^ = 0.515 and rji^^ = —2.809. Thus, 
we expect synchronization to occur at these values as rj is either increased or decreased away from 
zero. 



1.0n 




FIG. 4: Case E, different populations. Calculations of the order parameters versus rj confirm 
that the onset of synchronization occurs at ry^, = —2.809 and 0.515 (vertical lines), as predicted in 
Fig. ©. 
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FIG. 5: Case A, different populations. Panels as in Fig. Note that the standard branch cut 
leads to discontinutites and the occurrence of two roots for Im{r]^^^)(v) (solid lines, left panel), 
and none for Im{rj^'^^)(v) (dotted lines, left panel). From the right panel we find r/* = 2.189,4.501, 
taking care to obtain these from Re{ri^^^) (solid line, right panel). We expect to find synchronization 
onset at the smaller of these values. 



0.8n 




FIG. 6: Case A, different populations: The onset of synchronziation occurs at ry* = 2.189 (vertical 
line) , as predicted in Fig. ([5]) . No synchronization is observed for smaller values of rj (not shown) . 
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FIG. 9: (Color online) Three populations. The imaginary and real parts of r/*^^'^'^^ are plotted in 
blue, red, and green to illustrate the discontinuties due to branch cuts. The analysis proceeds as 
in the previous cases. Because of branch cuts, two roots occur on one on r/(^\ and none on 
r]^^\ Evaluating the real parts appropriately, we find r/* = —0.891, —0.564, 2.303. 




FIG. 10: Three populations. Synchronization occurs at rj^, = —0.564 and 2.303, as predicted in 

Fig. dHD. 



Tables 



Case 


Condition 




v* 


1 


T2 > AD 


-n 




2 


< AD 




4A 
T 


3 


D = 


-n 


2A 

T 


4 


T = 0, D <0 


-n 


_|_ 2A 


5 


T = 0, D>0 


no solution 


no solution 



16 



TABLE I: Solutions to Eq. ([8]) for two identical populations. D = det(K) and T = tr(K), where 
K is the connectivity matrix (Eq. ([I])), and A is the width parameter in Eq. ([6|). 



Case 


Matrix 


T 


D 




A 








1 


1 


4 


B 


( 




1 


-1 


1 


-4 


C 


( 


3 1 
-3.5 -1 


) 


2 


0.5 


2(2 - \/2) = 1.172 


D 


( 


^ -3 1^ 
V-3.5 IJ 


1 


-2 


0.5 


2(-2 + ^/2) = -1.172 


E 


( 


(-1 -1^ 
V 1 2 y 


1 


1 


-1 


-(l±^/5) = -3.236,1.236 


F 


( 


^ 1 1 ^ 

V-1 -2j 


1 


-1 


-1 


l±^/5 = -1.236,3.236 


G 


( 


^ 2 1 ^ 

V-3 -2 J 


1 





-1 


±2 


H 




(::;) 







1 


None 



TABLE IL Connectivity matricies K chosen to sample T — D space. 
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